Efficiency for preforming molecules from mixtures of light Fermi and heavy Bose 
atoms in optical lattices: the strong-coupling-expansion method 
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We discuss the application of a strong-coupling expansion (perturbation theory in the hopping) 
for studying light-Fermi-heavy-Bose (like 40 K- 87 Rb) mixtures in optical lattices. We use the strong- 
coupling method to evaluate the efficiency for pre-forming molecules, the entropy per particle and the 
thermal fluctuations. We show that within the strong interaction regime (and at high temperature), 
the strong-coupling expansion is an economical way to study this problem. In some cases, it remains 
valid even down to low temperatures. Because the computational effort is minimal, the strong- 
coupling approach allows us to work with much larger system sizes, where boundary effects can 
be eliminated, which is particularly important at higher temperatures. Since the strong-coupling 
approach is so efficient and accurate, it allows one to rapidly scan through parameter space in order 
to optimize the pre-forming of molecules on a lattice (by choosing the lattice depth and interspecies 
attraction). Based on the strong-coupling calculations, we test the thermometry scheme based on 
the fluctuation-dissipation theorem and find the scheme gives accurate temperature estimation even 
at very low temperature. We believe this approach and the calculation results will be useful in the 
design of the next generation of experiments, and will hopefully lead to the ability to form dipolar 
matter in the quantum degenerate regime. 



I. INTRODUCTION 

In recent years, there has been much interest in ultra- 
cold polar molecules as they have the promise for 
being a new state of quantum degenerate matter, with 
unique properties. In order to have a large dipole mo- 
ment, the polar molecules must be in their rovibrational 
ground state, where further cooling can ultimately lead 
to quantum degenerate dipolar matter Such polar 
molecules can have long-range, anisotropic or three-body 
interactions [3j , which may lead to novel quantum phases 
[3, [1[ and new applications in quantum information sci- 
ence @. In most ultra-cold polar molecule experiments, 
one starts with a mixture of ultra-cold gases of atoms of 
different species, for example various isotopic combina- 
tions of K and Rb These atoms can form a weakly 
bound state through a magnetic field sweep over the Fes- 
hbach resonance (ill Il2| . To create molecules with signif- 
icantly higher dipolar moments, the loosely bound Fes- 
hbach molecules are coherently transferred to a ground 
state with very high efficiency thr oug h stimulated Raman 
adiabatic passage (STIRAP) pHiq. 

Although the rate of transferring a Feshbach molecule 
to the ground state is very high, the overall efficiency for 
forming dipolar molecules is still low due to the low effi- 
ciency of forming the loosely bound Feshbach molecules 
during the field sweep. In Ref. [l^|, the fermionic 40 K and 
the bosonic 87 Rb atoms are trapped by an optical trap. 
The efficiency to form the Feshbach molecule depends on 
the phase-space density of the two species. But, because 
the Fermi cloud stops shrinking once it reaches the quan- 
tum degenerate regime, and the Bose cloud continues to 
shrink as it Bose condenses, this phase space density is 



low at low temperature, and never reaches appreciable 
sizes at higher temperatures, as the clouds become more 
diffuse. On the other hand, if the mixture is first loaded 
onto an optical lattice, the motion of the atoms can be 
more strongly confined, and it is possible to create a large 
area where exactly one atom of each species sits at the 
same lattice site, leading to a reduced three body loss 
[T3 | and almost unit efficiency [l7| for pre-forming the 
molecules. 

When mixtures of 40 K and 87 Rb are loaded into an op- 
tical lattice, the atoms of each species are influenced by 
the optical lattice differently [18j. With the same optical 
lattice depth, the heavy atoms usually have much lower 
tunneling rate than the light atoms because of their sig- 
nificantly larger mass. In Ref. [17], it was shown that 
for sufficient lattice depths, the hopping rate of Rb is 
more than an order of magnitude less than that of K. It 
is therefore reasonable to ignore the quantum effects of 
the tunneling of the heavy bosonic atoms while allowing 
the light fermionic atoms to hop between nearest neigh- 
bors (a classical effect of the motion of the Rb atoms is 
taken into account by averaging over all energetically fa- 
vorable distributions of Rb atoms) . Such systems can be 
described by the Fermi-Bose Falicov-Kimball model (l9l - 
[2l| . Using this model, we quantitatively determine the 
probability of having exactly one atom of each species per 
lattice site in order to optimize the formation of dipolar 
molecules. 

For the Falicov-Kimball model, the phenomena of pre- 
forming molecules has been discussed for Fermi-Fermi 
mixtures or Fermi-hard-core-Bose mixtures (22| on a ho- 
mogeneous lattice and Fermi-Fermi mixtures in a har- 
monic trap [23]. In previous work [l7| . we considered the 
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Fermi-soft-core-Bose mixtures in a harmonic trap and de- 
termined the efficiency for pre-forming molecules as the 
probability to have exactly one atom of each species per 
site. We used inhomogeneous dynamical mean-field the- 
ory (IDMFT) and Monte Carlo (MC) techniques to cal- 
culate the efficiency as well as the density profile and 
the entropy per particle. Both of these methods have 
advantages and disadvantages. The IDMFT approach 
is approximate for two-dimensional systems, but it can 
calculate both the efficiency and the entropy per parti- 
cle. The MC method is numerically exact after it reaches 
thermal equilibrium, but it can not calculate the contri- 
butions to entropy coming from the heavy particles. Both 
methods require large computational times to calculate 
properties of a trapped system of reasonable size. Using 
these methods, we have shown that the efficiency is sig- 
nificantly increased by first loading onto an optical lattice 
before forming the molecules and near unit efficiency can 
be achieved with parameters that are realistic for current 
experiments. 

The efficiency of pre-formed molecules is also likely to 
be affected by the heating (the temperature increase) in- 
duced by loading onto an optical lattice f24T-|26l| . Con- 
sidering that thermal fluctuations generally destroy the 
ordering and the localization of the particles, it is rea- 
sonable to expect that the efficiency of having exactly 
one Rb atom and one K atom per site should be re- 
duced if the temperature becomes too high. On the other 
hand, if the temperature is low enough, the presence of 
the lattice significantly increases the efficiency, almost 
to unity in the case of deep lattices. The temperature 
of the lattice system, however, remains difficult to mea- 
sure in experiment [27H3I] . Instead, it is often assumed 
that the process of loading atoms onto optical lattices 
is adiabatic and therefore the total entropy of the sys- 
tem is conserved [2414261 |32j, |33| . Determined based on 
the thermal properties of the gas before adding lattices, 
the entropy per particle is then used as an effective tem- 
perature scale for the lattice system (2f| |34| . There are 
also several proposals for directly determining the tem- 
perature for systems of bosonic atom s I35l -l37| , fermionic 
atoms [HI or the magnetic systems In Ref. (3l| . 

a general thermometry scheme is derived based on the 
fluctuation-dissipation theorem. Through quantum MC 
simulation, this proposal is shown to be applicable to 
the non-interacting fermionic systems (3f| and interact- 
ing bosonic systems (40| . 

In our paper, we discuss light-Fermi-heavy-Bose mix- 
tures in optical lattices based on the strong-coupling (SC) 
expansion method (perturbation theory in the hopping). 
The calculation is oriented to develop an efficient way 
of estimating the efficiency of pre-forming molecules for 
a given experimental system. With the SC expansion 
method, we obtain analytical expressions for the effi- 
ciency of pre-forming molecules, the entropy per parti- 
cle and the local charge compressibilities. The behavior 
of the efficiency is studied both as a function of entropy 
per particle and temperature. The determination of tem- 



perature is further studied by applying the thermometry 
proposal in Ref. [3l| to the Fermi-Bose mixture. To 
benchmark the accuracy, we compare the SC calculation 
with the IDMFT and MC calculations for all parame- 
ters considered. Overall, we find excellent agreement be- 
tween the three methods. Such agreement even extends 
to the low temperature region when the interaction is 
strong enough. This is particularly useful, given the fact 
that the SC expansion calculation is significantly faster 
than the IDMFT and MC calculations. Such a speedup 
makes it possible to consider much larger lattice sizes to 
eliminate the boundary effects, to scan the large param- 
eter space for optimal parameter regions for pre-formed 
molecules and to estimate the density fluctuations and 
other properties. 

The paper is organized as follows: in Sec. II, we dis- 
cuss the Fermi-Bose Falicov-Kimball model and define 
the efficiency for pre-forming molecules; in Sec. Ill, we 
discuss the formalism for evaluating the efficiency, the 
entropy, and other related quantities; in Sec. IV, we dis- 
cuss our result for various parameters and benchmark 
the SC expansion calculation with the IDMFT and MC 
calculations; in Sec. V, we discuss the application of the 
fluctuation-dissipation theorem for determining the tem- 
perature and we present our conclusions in Sec. VI. 



II. FERMI-BOSE FALICOV-KIMBALL MODEL 

For mixtures of heavy bosons and light fermions, such 
as 87 Rb/ 40 K mixtures, the hopping parameter for the 
heavy bosons ( 87 Rb) is usually more than an order of 
magnitude less than the hopping parameter for the light 
fermions ( 40 K) when one takes reasonable lattice depths 
(greater than 15 Rb recoil energies) [17j . In this case, we 
can ignore the quantum- mechanical effects of the hopping 
of the heavy bosons and describe such mixtures with the 
Fermi-Bose Falicov Kimball model in the presence of a 
trap potential. The Hamiltonian of this model is written 
as 

H = H a + H h = J2 Ho, + H h , (I) 

3 

with 

H 0j = (Vj-^fjfj + Uvfjfrbfa 

+ (V j -to)b]b j + ~U»bfa(b]b j -l), (2) 

and 

^ = -EV/]/i'' ( 3 ) 

Here, j, j' label the sites of a two-dimensional square 
lattice, with a lattice constant, a. The symbols /j and 
fj denote the creation and annihilation operators for the 
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fermions at lattice site j, respectively. The symbols fet and 
and bj denote the creation and annihilation operators for 
the bosons at lattice site j, respectively. The fermionic 
operators satisfy the canonical anticommutation relation 
{/ji/J'} = ail( i the bosonic operators satisfy the 
canonical commutation relation [6j , 6t, ] = djji. The 
quantity Vj is the trap potential, which is assumed to 
be a simple harmonic-oscillator potential centered at the 
center of the finite lattice. We assume that the jth site 
has the coordinate (xj ,yj), so that Vj can be written as 



V, = t 



Ml 
2to 



1 2 



(4) 



where f2 is the trap frequency. The quantity [if is the 
chemical potential for fermions and fib is the chemical 
potential for bosons. Combining the trap potential and 
the chemical potentials, we can define an effective posi- 
tion dependent local chemical potential for the fermions, 
fifj = fif — Vj, and for the bosons, /Lt^j = /ij — Vj. Ubf 
is the interaction energy between fermions and bosons 
and Ubb is the interaction energy between the soft-core 
bosons. The symbol —tjji is the hopping energy for 
fermions to hop from site j' to site j. We consider a 
general tjji for the formal developments in the earlier 
part of the next section, but later specialize to the case 
of nearest-neighbor hopping with amplitude t, which we 
will take to be the energy unit. We also set the lattice 
constant, a equal to one. 

The efficiency for pre-forming molecules is defined as 
the averaged joint probability of having exactly one boson 
and exactly one fermion on a lattice site, 



(5) 



with f3 = (ksT) 1 the inverse temperature. We define 
the projection operator P{ x for having exactly one boson 
and one fermion at site j, 

H,i = \ n b.j = hnf,j = = l,nf t j = 1|, (6) 

and N is the smaller value in the total numbers of bosons 
and fermions, Nf, and Nf. In our case, we assume equal 
number of bosons and fermions, therefore N = jVj = Nf. 

In general, one can obtain £ directly from Eq. ([5]) 
for a readily diagonalized Hamiltonian. In our case, the 
efficiency £ is derived by distinguishing the contribution 
from terms corresponding to rib : j = 1 in the expression 
for the density of fermions. We assume that the density 
of bosons and fermions at site j, pb^ and Pfj, can both 
be written as a series in terms of the occupation number 
of bosons at site j in the following way, 



Ph., 



(b]bj) =J2Wj(n b ,j)n b ,j, 



(7) 



Pf,j = (fjfj) = l>l W 3( n b,j) n fA n b,j)- 



(8) 



Here nb,j is the occupation number of bosons on site j, 
n b,j = Oj 1j •••• The coefficient V^jiribj) is the probability 
of having exactly nbj bosons at site j and the coefficient 
n/j (n{,j ) is the probability for having one fermion on site 
j for the occupation number n-bj. The joint probability 
of having exactly one boson and one fermion at site j can 
be written as 

£ i = w j( n b.3 = l )nf,j( n b,i = !): ( 9 ) 
and the efficiency £ is the average of £j over all sites, 



£ = 



N 



N 



(10) 



It can be shown that the expression for the efficiency ob- 
tained in this way is the same as from Eq. ([5]). Now, the 
efficiency is obtained directly from the density of bosons 
and fermions, which can be easily derived from the parti- 
tion function Z by taking derivatives with respect to the 
appropriate chemical potentials, 



Pb, 



1 ain(Z) 

P dfj,b,j 



and 



1 51n(Z) 



P dp/,. 



(11) 



(12) 



To study the behavior of the efficiency as a function 
of the entropy per particle, we evaluate the entropy per 
particle by dividing the total entropy by the total number 
of particles, 

s = S/(N b + Nf) 

k B hx(Z)-0k B ^p^/(N b + N f ). (13) 

It is worthwhile noticing that the formalism develop- 
ment in this section is based on the grand-canonical en- 
semble. This ensemble is appropriate because we assume 
that in the lattice system both the energy and the num- 
ber of particles fluctuate. This may seem in contradiction 
with the use of the entropy per particle as an effective 
temperature scale, because strictly speaking entropy is 
used as a parameter only for the micro-canonical ensem- 
ble. This contradiction is resolved because the entropy 
per particle is assumed to be conserved during the pro- 
cess of turning on the optical lattice. It is a conserved 
quantity when comparing the systems before and after 
turning on the optical lattice, which is particularly use- 
ful from the experimental point of view, since the exper- 
iments often start without the optical lattices. For the 
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lattice system itself, assuming it is in thermal equilib- 
rium, it is more reasonable to consider it with the grand- 
canonical ensemble, allowing the energy and number fluc- 
tuations. The difference between the different ensembles 
of course is not problematic if we assume the system is 
large enough to be in the thermodynamical limit, where 
all three ensembles are equivalent. 



III. STRONG-COUPLING EXPANSION 
FORMALISM 



In this section, we explain the SC expansion formalism. 
We first discuss the evaluation of the partition function 
Z, approximated by the second-order expansion in terms 
of the hopping, Hh- From the partition function, we de- 
rive the expressions for the density of fermions shown in 
Eqs. (1211) to (JMJ), the density of bosons in Eqs. (J35J to 
(1571) . the efficiency in Eqs. d35]) to (gD]) and the entropy 
per particle in Eqs. (|4"4")l to (14"7) . For readers who pre- 
fer to see the final expressions, we suggest skipping the 
following derivation and referring to the equations listed 
above for the corresponding quantities. 

The evaluation of the partition function in the SC 
approach starts with the exact solution of the atomic 
Hamiltonian Hq. Hence, we use an interaction picture 
with respect to Hq, where for any operator .4, we de- 
fine the (imaginary) time-dependent operator A{t) — 
e rHo Ae~ rH ° . The partition function is written using the 
standard relation, 



Z = Tr (e~ pB ) = Tr ( e -^ H °U(i3, 0)) . 



(14) 



Here, U(/3, 0) = T r exp J Q Hh(r)dT is the evolution op- 
erator with T T being the time-ordering operator for imag- 
inary times. Expanding the exponential in U(/3,0) up to 
second order in Hh (t) and evaluating the resulting traces 
with respect to equilibrium ensembles of Hq, we have 

U((3,0) ~ 1 + 



dr x \ dT 2 T T H h (n)H h {T 2 ). (15) 



Here, the first order correction to the partition function 
vanishes because the hopping connects different sites. 
Substituting Eq. (fTSJ) into Eq. (TH)) , we can write the 
partition function as, 



z = z(°)(i + z< 2 )), 



(16) 



where Z^ is the partition function in the atomic limit 
(t = 0), 

= Tr (e- tiHo ) , (17) 
and Z^ corresponds to the second-order term in the 



expansion of U divided by Z^°\ 



Z^ = 



2Z(°) 



Tr 



-PH 



/ / dnd^TrHhiT^H^] 

Jo Jo 



(18) 

To simplify the notation, we introduce fif,j{nb,j) to 
represent the negative of the fermionic part of the Hamil- 
tonian Hqj [Eq. ([2])] when there is a fcrmion at site j, 



/</,./('"-.; 1 /'; Vj - UbfUbj, 



(19) 



and Hb,j{nb,j) for the negative of the bosonic part of the 
Hamiltonian Hqj, 

= (pb - Vj)nb,j - U bb n btj (nb t j - l)/2. (20) 

Both ft, f j and fib ,j depend on the number of bosons at 
site j. The effective fugacities for bosonic and fermionic 
particles can then be written as the exponential of ftfj 
and p b ,j respectively, 



and 



i>f,j(m) = ex P \fi^f,j(n b ,j)] 



<Pb,j(n b ) = exp \PHb,j(nb,j)] ■ 



(21) 



(22) 



The atomic partition function Z^ can then be written 
in terms of the effective fugacities as, 



(23) 



where Z^ is the atomic partition function at site j, 



(24) 



Now we evaluate the second term in the partition func- 
tion, Z^ of Eq. (j3"Tj) . To satisfy the total number 
conservation, only terms with j = kl and j' = k in 
Hh(T\)Hh{T2) are non-zero after the trace and Z^ is 
reduced into a sum of products of the fermionic annihi- 
lation and creation operators at the same site, 



1 f p f p 
Z<2) = o / / dndT 2 y2t jk t k 
IJo Jo jk 



xTr 
xTr 



AO) 



Tre-^fUn)/^)] /Z 



T T e-^f k (n)fl(T 2 ) !Zf. (25) 



Using the cyclic permutation relationship of the trace, 
the products can be represented by the local atomic 
Green's function, 



Gn(r) 



-Tr 



T T e-^f 3 {T)f}{Q) /Zf\ (26) 
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and Z^ 2 ' is expressed as integrations of the atomic Substituting Eq. (I2TJ1) into Eg. (|26p . we obtain the atomic 
Green's functions in terms of their relative times, Green's function in terms of the effective fugacities as, 

,/3 ,/3 



1 f p f 

-2 (2) = -- / / dTidr 2 tjktkj G k k (n -T 2 )Gjj (t 2 -ti). 

(27) 

Solving the Heisenberg equation of motion for the an- 
nihilation operator fj(r), 



djjir) 
dr 



= e HoT [HoJ^e- 11 ^ 



(28) 



Gjj(r) = 



V <t>b,j{nb) mi j > n 

j 

<j>b,i( n b)4>S,i( n b) 



<fb,jY"b)Vf,jY"b) TUf ,■ _ ^ A 



(30) 



one easily finds the expression for the annihilation oper- 
ator /j(t) in the interaction picture, 



/i(r) = e^^)V J -(0), 



(29) 



We now perform the integration over n and T2 in Z^ 
and obtain the final expression for Z^ 2 \ 



n 6, j ™b,fc 



4>b,]{n b ^)4>b,k{ n b,k) P [0/,j("b,j) - 0/,fc(»b,fc)] 

Z (0) Z W ftf,j(n b ,j) - fif,k( n b,j) 



(31) 



Note that the partition function we derived here is not 
limited to the case of nearest-neighbor hopping with a 
uniform hopping parameter t. Eq. (|31[1 can be applied 
to describe hopping between arbitrary sites j and k and 
the hopping parameter tjk can vary over different sites of 
the lattice. 



Observables are evaluated by taking appropriate 
derivatives of the partition function. In calculating the 
derivatives, we truncate all final expressions to include 
only terms through the order of t 2 k . Also note that be- 
cause sites j and k are different sites, we do not nor- 
mally have denominators equal to zero in Eq. (|31j) . but 
in any case, the formulas are always finite as can be veri- 
fied by l'Hopital's rule. During numerical calculations of 
the observables, the denominator, p-f,j(nb,j) — P>f,k(nb,j), 
can become too small and cause numerical errors. In 
our calculations, we use the Taylor expansion in terms 
of pf,j(nb t j) — Pf,k{ n b,j) around zero when the absolute 
value of jlf,ji. n b,j) — fif,k(nb,j) is l ess than 10~ 5 . 



where pf j is the density of fermions in the atomic limit 



p {0) - 

Pf 



1 p dfif 



(2) 

and ■ is the total contribution to the density at site j 
from particles hopping from all possible sites, 



r(0) 



(33) 



(2) 



^t 3k tkj \ 



0bj(«-&j)0b,fc(«-b,fc) 



o.(0) _(0) 
3 k 



P 



C 1 - pfl)<t>fA n h,j) + Pjl<t>f,k{nb,k) 



(/>f,k( n b,j) ~ <l>f,j(nb,k) 
iPf,j(nb,j) - P-fA n b,j)Y 



(34) 



Similarly, the density of bosons at site j is written as a 
sum of the atomic density and the hopping contribution 

as, 



The density distribution is evaluated by taking the 
derivative of the partition function with respect to the 
appropriate local chemical potential [Eqs. (fT2| and (|lip]. 
For the density of fermions at site j, the expression con- 
stitutes two terms corresponding to derivatives from Z^> 
andZ( 2 ). 



Pf,3 



1 dln(Z) 

P d Pf,3 



(0) 
Pf'j 



(2) 
Pf,P 



(32) 



Pb,j 



where 



"> J P dm- 



co: 



1 dln(Z) 

P C>Pb,j 



o (0) + P (2) 

Pb,j ^ Pb.j ' 



AO) 



(35) 



E„ b , n b,j<f>b,j{rib,j) [1 + <f>f,j(nb,j)] 



(36) 
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and 



.-,(2) 



/'/...; = ^tjktkj ^ 

P[4>f,j(nb,j) ~ 4>fA n b,k)] 



, (Q)AbA n bA$bA n b,k) 

[nb,j - Pbj)— 



z (o) z (o) 

3 k 



(37) 



The expression for the efficiency is obtained from the 
density distributions of the fermions and bosons. Similar 
to the expression for the densities, the efficiency consists 
of two terms, one corresponding to the atomic limit and 
one corresponding to the contributions from the hopping, 



& = £ 



(o) 



£ 



(2) 



where 



and 



g (0) = (j>b,j(n b , j)(l>f,j(nb,j)\, 



AO) 



(38) 



(39) 



7(2) 



4>b,j{nb,j)4>b,k(n b , k ) ^ 



z (o) z (o) 

j k 



^b.j(n' btj )(l)fj(n' bj ) 



AO) 



.P[<t>fA n b,j) - 0f,k{ n b,k)} 



+5. 



"6,1,1 



WjKj) - fJ-f,k( n b,j) 

</>/,fc(»bj) ~ (j>f,j{n b ,k) 

\fifA n bj) - vfA n b,j)Y 



(40) 



For the trapped system, the local chemical potential 
]Xj includes both the global chemical potential /i and the 
trapping potential Vj . The derivatives with regard to the 
local chemical potential or the chemical potential leads 
to different physical quantities. For the Fermi-Bose mix- 
ture considered here, the cross-derivatives should also be 
evaluated. Specifically, the derivative with regard to the 
global chemical potential (fi b + fif) corresponds to the 
total number fluctuations, 



d 2 hxZ 



P 



((N f + N b ) )-(N f + N b y 



(41) 



Here we define the total number operators, Nf 



_tifi and N b = Y^A 
ity is introduced as the response of the local density to 



J2j f] fj an d N b — J^jbjbj. The global compressibil- 



the change of the global chemical potentials, 

d 2 lnZ 



d(fJ,f,j + M6,j)<9(>6 + M/) 



■if fi+b}b 3 )(N f +N b ) 



(42) 



And the local compressibility, or the onsite number fluc- 
tuation, is determined from the derivatives with regard 
to the local chemical potential, 



d 2 \nZ 



' yJ U'b.j ■ iu.j) 



P 



(U Vi + bfo) )-(f]fs + WY 



(43) 



Both the global and local compressibilities are derivatives 
of the density distributions and can be obtained from the 
density expressions above. 

Finally, we obtain the expression for the entropy per 
particle defined in Eq. (fTU|) by averaging the total en- 
tropy of the system and we again write the entropy per 
particle in terms of the atomic limit expression and the 
contributions from the hopping, 



j 



(2) 



(44) 



Here Sj is the entropy at site j in the atomic limit, 



AO) 



(45) 



where the parameter e-j corresponds to the onsite energy 
at site j in the atomic limit, 

<91n(Z(°)) 
€] ~ dp 



75j E {^,j( n b,j)^b,j(n b ,j) [1 + 4>fA n bA\ 
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3 n b,J 



(46) 



The averaged contributions from the hopping at site j is 

S, 



(2) 



<i>b,j(n b A<l> b ,k( n b,k) 



= --E E -o ,. 

x {[feKj) ~ <Pf,k(n b , k )} 

x [p- b ,j(n b ,j) + fL b ,k(n b ,k) - tj - £fe] 

+ P>fA n bA<t>fA n bA - M/^K,fc)</ > /,fe(«bj)}(47) 

This ends the discussion on the derivation of the SC 
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expansion method formulas. In general, the expressions 
obtained above are accurate in the case when the hop- 
ping is much smaller than interaction strength and the 
temperature is very high (J3t is small) . In this parameter 
region, the SC method can evaluate physical quantities, 
like the density distribution, efficiency, compressibility 
and entropy, very efficiently. The total number of parti- 
cles is fixed by varying the chemical potentials, //& and 
(if. To maximize the efficiency and reduce three body 
loss, we consider the low density region with attractive 
interspecies interactions and repulsive bosonic interac- 
tions. For other strong-coupling regions, the formulas 
developed above are equally applicable but not further 
discussed in this paper. 



IV. RESULTS 

A. Comparison with the IDMFT and MC 
calculations 

For a perturbative method like the SC expansion 
method, it is always necessary to determine the param- 
eter regions where the approximation is valid. Here, we 
use the previous results obtained from IDMFT and MC 
methods [It} as a reference to determine the accuracy of 
the SC calculation. It is also worthwhile to notice that 
the three methods require substantially different compu- 
tational times. The SC calculation usually takes less than 
1 CPU hour while for the same system the IDMFT calcu- 
lation takes on the order of 10 5 CPU hours. We consider 
all the parameters used in the previous work [TJ. The 
lattice is 50 x 50 square lattice with the trap frequency O 
for both species fixed at Ml / 2ta = 1/11. The parameters 
Ubf and Ubb are chosen based on a typical experimen- 
tal setup : Ubf It = -8, -12, -16 for U bb /t = 11.5 and 
Ubf It = -2, -6, -10 for Ubb = 5.7. The total number of 
bosons and fermions are set to be 625. We consider the 
temperature range 0.05i/&B to 20i/fcg. 

In Fig. [TJ we show the efficiency as a function tem- 
perature calculated with the three methods. Overall, we 
find excellent agreement between the SC result and the 
IDMFT and MC calculations and it is clear that high 
(unit) efficiency can be achieved when the temperature is 
low (T ~ O.lt/kg) and the interaction is large compared 
with t. In the case of Ubf = —2t and Ubb — 5.7t, the 
SC calculation starts to deviate greatly from the IDMFT 
and MC calculation when T < lt/ks- It is worth noting 
that for T > lt/ks, the SC calculations agree nicely with 
the other methods even for a relatively weak cross-species 
interaction, Ubf — — 2t. 

The difference between the SC calculation and the 
other two methods can be understood from the fact that 
the SC method is a perturbative method based on the 
atomic limit of the Hamiltonian, t = and that the prop- 
erties derived from the SC expansion are dominated by 
the atomic-limit behavior with relatively small correc- 
tions from the hopping. In the atomic limit, bosons and 



fermions are completely localized and the only density 
fluctuations are due to thermal fluctuations. For the low 
density case considered here, the bosons always form a 
plateau of unit filling at the center of trap at low tem- 
perature and the fermions are attracted by the bosons 
one by one and form an almost identical plateau. The 
efficiency therefore always converges to unity as temper- 
ature deceases. In Fig. [TJ we indeed find the efficiency 
from the SC calculation always goes to one at low tem- 
peratures. The convergence to unity is also true for the 
IDMFT and MC calculations for all the parameters ex- 
cept for U\,f It — —2 and Ubb/t = 5.7. That's where the 
SC calculation differs from the IDMFT and MC calcula- 
tion. It is reasonable to assume that the SC calculation 
can be applied to the region where the ground state of 
the system is a localized, Mott insulator like state. 

The SC calculation of the entropy per particle is also 
compared with the IDMFT and MC calculations for all 
the parameters using Eqs. (|44M47[) . The conclusion of 
the comparison is similar with the efficiency calculation, 
that the SC calculation is accurate except for U b f = —It. 
In Fig. [5] we use one example, Ubb = H-5t and Ubf = 
— lQt, to represent all the cases where the SC calculation 
agrees with the IDMFT calculation. As the temperature 
increases, the entropy per particle starts to saturate at 
around ~ 2.3fcs. In the next section, we will show that 
this saturation is actually the result of finite-size effects. 

In Fig. [31 we show the behavior of the efficiency as a 
function of the entropy per particle. This figure can be 
compared with Fig. 2 in Ref. [13, where the IDMFT 
calculation is discussed. We verify the findings from the 
previous work that for strongly attractive inter-species 
interactions, an efficiency of 100% can be achieved at 
low temperature (low entropy) region. For an entropy 
per particle around 1/cb, a 80% efficiency can still be 
reached. This efficiency is much higher than what has 
been achieved in experiment [l5j . 

In the following discussion on the SC calculation re- 
sult, we no longer consider the case of Ubf — —It. This 
is also based on the consideration that the interaction of 
Ubf = — 2t is too weak to achieve the desired high effi- 
ciency of pre-formed molecules and therefore is not in the 
parameter region of the main interest in this paper. 



B. Finite-size effects 

In our calculations, we always assume a hard-wall 
boundary condition at the edge of the lattice. In ex- 
periments, however, the atoms are confined only by the 
trapping potential. This additional confinement imposed 
by the boundary condition can potentially affect the ac- 
curacy of our calculation. This finite-size effect can be 
neglected if the system is so large that the atoms trapped 
by the trapping potential almost never reach the edge 
the system. This, however, is not always the case for the 
50 x 50 lattice. This problem is difficult to address with 
the IDMFT and MC methods, because of the high com- 




Figure 1: (Color on-line) Efficiency £ as a function of temperature calculated by the SC (red cross), IDMFT (blue triangle) and 
MC (green circle) methods. The interaction parameters, Ubb and Ubf, are shown in each plot. In (a), the SC calculation differs 
from the other two methods for T < lt/ks- For this region, the SC expansion formulas derived here are no longer accurate. In 
(b)-(f), all three methods give almost identical results. These calculations also show that almost 100% efficiency is reached for 
relatively strong attraction, Ubf > — 6t, at low temperature, T < t/ks- 
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Figure 2: (Color on-line) (a) Entropy per particle as a func- 
tion of temperature T. The SC calculation is marked with 
red crosses and the IDMFT calculation by the blue line. We 
find excellent agreement between the SC calculation and the 
IDMFT calculation. 

putational costs. The SC method, on the other hand, 
can calculate much larger systems for a fraction of the 
cost. 

In this section, we discuss our calculation for different 
lattice sizes and discuss finite-size effects for different lat- 
tice sizes. To benchmark the SC calculations, the trap 
frequency and the total number of particles are fixed for 
all different lattice sizes. We assume the largest lattice 
sizes are sufficient to neglect the finite-size effects. In Fig. 
[2J we show the density profile as a function of the lattice 



size at two temperatures, T = It/ks (a) and T = 20t/ks 
(b). Here Fig. 0] (a) represents the scaling behavior in 
the low-temperature region, where there is no significant 
difference between different lattice sizes and Fig. 2] (b) 
represents the scaling behavior in the high-temperature 
region, where the system of small lattice size is highly 
affected by the boundary effect. Note that the horizontal 
axes are different scales in the two panels. The param- 
eters used in the plots are Ubf = — 16t and Ubb = 11.54. 
We find similar behavior of the density profile for all the 
other parameters. 

In Fig. [S] (a), we show entropy per particle as a func- 
tion of temperature at different lattice sizes. In this plot, 
we find that for small lattices, the entropy per particle 
becomes saturated at high temperature, while for large 
lattices it keep increasing as the temperature increases. 
The saturation is understood as the result of the finite- 
size effects. When the temperature is high, atoms tend 
to expand to a larger area in the trap, which leads to 
a large cloud size and higher entropy. When atoms ex- 
pand to the edge of the lattice, the possible occupied sites 
are now constrained and the entropy stays similar even 
though the temperature increases, hence the saturation. 
When the lattice is sufficiently large, atoms can freely ex- 
pand as the temperature increases and the entropy keeps 
increasing. 

The confinement of the atomic cloud in high tempera- 
ture also affects the efficiency calculation. In Fig. [5] (b), 
we find that the efficiency saturates to a higher value 



9 




U =-6t (IDFMT) 



0.5 1 1.5 2 2.5 
Entropy per particle (kj 

Figure 3: (Color on-line) Efficiency as a function of entropy 
per particle for different interaction parameters. Note here 
that we didn't include the case of Ubf = —24, because it is 
already shown in Fig. [T]that the SC calculation is not accurate 
for low temperatures in this case. In (a) and (b), we consider 
two different bosonic interaction strengths and five different 
inter-species interaction strengths. For all parameters, the 
efficiency reaches 100% when the entropy is very low. For 
an intermediate entropy, with an entropy per particle around 
lfcs, the efficiency is around 80%. 
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Figure 4: (Color on-line) Finite-size effects on the radial den- 
sity profile. We assume a two-dimensional N x N square 
lattice with hard-wall boundary conditions. The dotted lines 
indicate the boundaries of different lattices. The interaction 
parameters are Ubf = —164, Utb = 11.54. We use the density 
distribution of the fermionic particle to represent the general 
dependence of density on the lattice size. In (a), we consider 
the case of low temperature, T = t/ks- Here, the density dis- 
tribution is concentrated at the center of the trap and there 
is no difference between different lattice sizes. In (b), we con- 
sider the case of high temperature at T = 204/fcs. Here, 
the density is confined mainly by the size of the lattice. For 
N — 50, the density is confined at the edge of the lattice, 
r — 25. For N = 100, the density is again confined at the 
edge, r = 50. For both N = 200 and 300, the density goes to 
zero before reaching the edge of the lattice and the two distri- 
butions overlap with each other. We estimate that finite-size 
effects are eliminated for the 300 x 300 square lattice for the 
trap frequency and number of particles considered here. 



for smaller lattices. This is because the confinement in- 
creases the density overlap between the two species. In 
the low temperature region, the atoms are close to unit 
filling at the center of the trap and the efficiency is similar 
for all difference lattice sizes. 



V. THERMOMETRY 

A. Temperature and Density fluctuations 

Based on the fluctuation-dissipation theorem, the com- 
>ressibility can be related to the density fluctuations as 



We find that a lattice of 300 x 300 sites is sufficient 
to eliminate the finite-size effects for our parameter re- 
gions. Hence, we use this lattice size for the efficiency 
and entropy per particle calculations. In fig. [SJ we show 
the result for the efficiency as a function of the entropy 
per particle. We estimate the calculation result from the 
50 x 50 lattice is accurate when the temperature is around 
or below T = 1.254/fc B . 



« = = ^ i(p(r)N) - (p(r))(N)] , (48) 

where p(r) is the radial density profile, /x is the chemi- 
cal potential and N is the total number of particles. For 
a system with a spherically symmetric harmonic trap- 
ping potential, — V t r 2 , the local chemical potential at a 
radial distance r is fi — Vtr 2 . Within the local density 
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Figure 5: (Color on-line) Finite-size effects on the entropy 
per particle and the efficiency. We assume a two-dimensional 
NxN square lattice with hard- wall boundary conditions. The 
interaction parameters are Utf = — 16t, Ubb = 11. 5t. In (a), 
we show the behavior of the entropy per particle as a function 
of temperature for different system sizes. We see the entropy 
is significantly affected by the finite size when the lattice is 
smaller than around 200 x 200. The finite-size effect is not 
noticeable at lower temperature (T < It/fee). 




1 2 3 4 5 




1 2 3 4 5 
Entropy per particle (k ) 



Figure 6: (Color on-line) Efficiency as a function of entropy 
per particle for a 300 x 300 square lattice system. We consider 
625 atoms for each species. Compared with Fig. O the effi- 
ciency is significantly higher for the same value of the entropy 
per particle in the 300 x 300 lattice system when the entropy 
per particle is large. On the other hand, the behavior is sim- 
ilar in both lattice systems when the entropy per particle is 
less than Iks- The unit efficiency is reached roughly when 
the entropy per particle is less than 0.5fcs. 



approximation, the trapping potential is interpreted as a 
variance in the chemical potential and the compressibil- 
ity in the trapped system can be re- written as a function 
of the density gradient, 



dp{r) 
dp 



1 dp(r) 
'ZVtX dr 



(49) 



These two equations lead to a simple relationship be- 
tween the density gradient and the density fluctuations 
in the trapped system, 



1 



1 dp{r) 
2V^ dr = I^f 



[(p(r)N)-(p(r))(N)} 



(50) 



Based on this relationship, one can determine the tem- 
perature from the independently measured density gra- 
dient and density fluctuations. For a two dimensional 
system, a simplified relationship can be found by inte- 
grating the above equation over all the two dimensional 
plane, 



7T , . 1 

VP = TT^F 
V t k B i 



«iV 2 > - (TV) 2 



(51) 



Here, p(Q) stands for the density at the center of the trap. 



With the development of in situ measurements, it is 
now possible to measure the density gradient and the 
fluctuations [H, |4l| in experiment and this thermome- 
try scheme has shown promise to be a reliable way of 
estimating the temperature (3ll . |40| . Here we test this 
method for the Bose-Fermi mixtures and Eqs. (|50l and 
I51| are extended to mixtures by considering the density 
as the total density of both species and the total num- 
ber as the total number of both species. With the SC 
method, we calculate the density gradient directly from 
the density profile expressions. To simulate the fluctu- 
ations measured in the experiments, we use a simplified 
MC simulation explained in the next section. 
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Figure 7: (Color on-line) Density fluctuations averaged over different numbers of samples. The density fluctuations shown 
here are the total density fluctuations divided by the input temperature, T = 2tk^ j 1 . All the fluctuations are compared with 
— (2Vtr)~ 1 d(pb + p f) / dr . According to Eq. (|50|> . these two quantities should be equal to each other. In (a)-(c), the total number 
of configuration generated is 2 x 10 5 , with a different sampling strategy. In (a), one sample is taken at every 10 3 configurations, 
which gives a total of 200 samples to average over. The statistical error in this case is very large. In (b), one sample is taken 
at every 100 configurations, which gives a total of 2000 samples. The statistical error is reduced compared with (a). In (c), the 
total number of samples is 2 x 10 4 . The statistical error is the smallest among (a) to (c). In (d), a total of 2 x 10 6 configurations 
are generated and 2 x 10 4 samples are taken at every 100 configurations. 



B. Fluctuation calculation 



probability in the atomic limit, 



The MC simulation method generates a large collec- 
tion of states (or configurations) that satisfies the thermal 
equilibrium criteria. Such collection of states constitutes 
a thermal ensemble. In the ensemble, each state (or con- 
figuration) gives one density distribution, analogous to 
one single shot image of the density in the experiment. 
By averaging over all configurations, one obtains the av- 
eraged distribution of particles. Deviations between dif- 
ferent configurations are the fluctuations. In our sim- 
plified MC method, we use the SC method to determine 
the density distribution for a given temperature and then 
use the probability as a reference for configuration gen- 
eration. The ensemble of configurations is decided to be 
large enough if it can reproduce the input probabilities. 



pj-(0) _ 0b,j(n)(j)fj(n) 



and the contributions from the hopping, 



P, 



i(i) 



71,1 



4>b,j( n )<l > f,j( n ) 



k 

E 

„ - 0fA n b,k) 



4>b,j{nb,j)4>bA n b,k) 
o.(0)o.(0) 
j k 



M/jKj) - M/,fc(»6,fc) 



(53) 



(54) 



Determining the joint probability : the joint probabil- 
ity, P£ m , is the joint probability of having n bosons and 
m fermions at site j. For m = 1, the joint probability 
of having n bosons and 1 fermion at site j can be found 
from the fermionic density distribution, similar to the 
calculation of the local efficiency £j (indeed, £j — Pd), 



pj _ pi(o) _ pjy-/ 1 pj\^ 

r n.\ ~ r n.\ r n.\ ^ r n.\ 



J(l) 



>J(2) 



(52) 



where we again write the probability as a sum of the 



P. 



n.l 



y-v 4>b,j(n)4>b,k{ n b,k) x 



+ 



7 (0) o-(O) 
k n b ,k ^k 

4>fA n ) - 4>fA n b,k) 



\p>fA n ) - t*fA n b,k)Y 



P<t>fA n ) 



M/jfa) - VfA n b,k) 

(55) 



Once the joint probability f~ 1 is determined, the com- 
plementary probability P^ is found based on the rela- 
tionship in the atomic limit, 
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E 



pi(o) 



</>6,j(n) 



Z 



(0) 



= 1. 



(56) 



Taking into the account the hopping contributions, we 
can write P 3 n as, 



- r n,0 ~~ 



(0) r n,l 



J(2) 



(57) 



We assume each lattice site is independent and the joint 
probability at site j is sufficient to determine the density 
distribution at site j. The joint probabilities are evalu- 
ated for all the lattice sites and stored in a table before 
the MC procedure. 

Simulation procedure: we use a random number gener- 
ator to generate configurations with reference to the joint 
probability table. Specifically the simulation includes the 
following steps: 

1) Create a table for the values of P£ m corresponding 
to the sum of the joint probability of having up to n 
bosons and up to m fermions at site j = 1, i.e. 



P J , 



EE 

k=0 1=0 



P 



(58) 



2) Generate a random number x between and 1. 

3) Find the smallest Pi, , that is larger than x. The 
number of bosons and fermions at site j is then equal to 
n' and ml . 

4) Repeat steps (2) and (3) to another site, j — 2, until 
all the lattice sites are considered. Store the configura- 
tion. 

5) Repeat steps (2)-(4) M times to generate J\f config- 
urations. 

To avoid auto-correlation between adjacent configura- 
tions, we choose every other M 3> 1 configurations as 
samples. The total number of samples is then N s — 
N I M.. Averaging over all the samples, we obtain the 
fermionic and bosonic part of the density fluctuation as 

5 f(b){r) = (p m (r)(N f + N b )) - (p f{b) (r)) (N f + N b ), 

(59) 

and the total density fluctuation is the sum of 8f and S b . 
The total number fluctuation is defined as 



((N f + N b ) 2 ) - (N f + N b ) 2 



(60) 



Here the bracket stands for the averaging over all sam- 
ples in analogy to the experimental measurement of the 
fluctuations. 
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Figure 8: (Color on-line) Extracted temperature as a function 
of the input temperature. The fluctuations are obtained from 
2 x 10 4 samples out of 2 x 10 6 configurations. The value of 
Ti is the mean of Ti(r) averaged over 12 < r/a < 25 and 
the error-bar for Ti is the standard deviation in Ti(r) [Eq. 
(pTj) ]. The value of T 2 is obtained through Eq. (55). The 
input temperature T is drawn as a straight blue line in both 
plots. In (a), we show our result for the full range of the input 
temperature, from T = 0.2£/fcs to 20t/fc_e. In this plot, we 
find very good overall agreement of Ti and T2 with the input 
temperature for the temperature range considered, particu- 
larly for T > lt/ks- In (b), we blow-up the area inside the 
black square in (a), which corresponds to the low temperature 
region, where T = 0.2£/fcs to 0.5t/fcs. In this region, we find 
that Ti shows large relative fluctuations (deviation) from the 
mean value and the mean value of Ti differs relatively greater 
from T. The extracted temperature T2 however still shows 
excellent agreement with the input temperature. 



C. Results 

The fluctuation calculation is carried out for a 300 x 300 
lattice with all five sets of parameters. Overall we find 
very similar behavior for all the parameters and we use 
parameters U bb — 11. 5t and U b f = — 16< as an exam- 
ple. In our simulation, the fluctuations between different 
configurations are from both the random number gen- 
erator and the thermal fluctuations. The difference be- 
tween them is that the thermal fluctuations are indepen- 
dent of ensemble sizes and the sampling size. We find 
that the correct thermal fluctuation calculation requires 
a large number of samples (~ 10 4 ) and large ensemble 
sizes (~ 10 6 ). Because of the similarity between the sim- 
ulation and experimental measurement, this may suggest 
that a large number of shots are needed in the experi- 
ments to obtain the correct thermal fluctuations. Note 
that we consider here the results from a single plane as 
one shot, not the averaged results over many planes as 
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reported in Ref. (io| . 

In Fig. [JJ we discuss the sampling effects by com- 
paring the fluctuations obtained from different samples 
with the compressibility calculated from the density gra- 
dient (2Vtr)^ 1 d [p/(r) + Pb(r)] /dr. When the number of 
samples are small, the fluctuations are largely random 
deviations from the average value. In Fig. [JJ (a), the 
fluctuations can be equally positive and negative, which 
does not even satisfy the condition that the total fluctu- 
ations should be always positive. As the number of sam- 
ples grows, the random noise starts to be averaged out 
and the fluctuations start to agree with the fluctuation- 
dissipation theorem. In Fig. [JJ (c) , the fluctuations agree 
very nicely with the relationship predicted by Eq. ([50]) . 
To show that 10 4 samples are sufficient, we consider an 
even larger ensemble, with 2 x 10 6 configurations [Fig. 
[7j(d)] and find that the two ensembles produce almost 
identical results. This shows that the fluctuation cal- 
culations obtained in this way are independent of the 
ensemble size and should correspond to the thermal fluc- 
tuations of the system. 

With Eqs. ([50)) and (p)Tj) , we define two extracted tem- 
peratures, T\ and T%. Let T\(r) be the extracted temper- 
ature obtained in terms of the density fluctuations and 
the density gradient, 



k B T{(r) 



S f (r)+5 b (r) 



{2V t r)- 1 d(p f (r) + p b (r))/dr 



12 < r/a < 25. 

(61) 

Here, we choose the radial distance to be larger than 
12 lattice sites because the quantity (2V t r)^ 1 d(pf (r) + 
Pb(r))/dr diverges as r — > for a finite density gradient 
and for small r, it goes to zero as the density develops 
a plateau at unit filling at low temperature. The radial 
distance is less than 25 lattice sites, because the density 
is almost zero in the outer regions and that increases the 
relative error. Together, we find that r between 12 and 
25 sites to be the best region to fit the fluctuation and 
the compressibility with each other. We also note that 
Eq. ([50)) still holds if one considers just the fermionic or 
bosonic part of the system, i.e. keep only the index / or 
b in S and p. 

The temperature T 2 is obtained based on Eq. (1511) . 
which translates into the following expression for the 
Bose-Fermi mixture, 



We also find in the low temperature region, T2 generally 
fits better with the input temperature [Fig. [5] (b)[. This 
finding suggests that the statistical error introduced by 
the numerically generated ensemble is lower in the cal- 
culation of T2 and this could be because the calculation 
of T2 only involves the first-order observable, the den- 
sity and the total number fluctuation, whereas, for the 
calculation of T\, we need to calculate the second-order 
observable, the density fluctuation, which may be more 
susceptible to statistical errors in the numerical simula- 
tion. 



VI. CONCLUSION 

The SC expansion method is a very efficient way of 
studying thermal properties of strongly interacting sys- 
tems. Through comparison with the IDMFT and MC 
calculations, we show that the strong coupling expansion 
method can be used for a wide range of parameters even 
at low temperature when the attractive interaction be- 
tween the two species is relatively strong. We use the SC 
method to evaluate the finite-size effects in our previous 
calculations. This leads to important modifications of the 
efficiency and the entropy per particle at high tempera- 
ture. The SC calculation also provides a way to simulate 
experimental measurements of the fluctuations. Based 
on the simulation, we find that the thermometry pro- 
posal based on the fluctuation-dissipation theorem is ac- 
curate for heavy-bose-light-fermi mixtures. This scheme 
suggests an effective thermometry scheme that works in 
the extreme low temperature in the deep lattice region. 
Overall, our work shows a promising way of creating 
strongly interacting quantum degenerate dipolar matter 
by loading the mixtures onto an optical lattice before 
the molecule formation. In addition to higher efficiency, 
the molecules created in this way are already situated in 
the optical lattice and can be directly adjusted to realize 
the novel quantum phases that require the presence of 
a lattice. It is also worth noting that the SC approach 
can be used to study other mixtures with modifications. 
For Fermi-Fermi mixtures like 6 Li- 40 K, it would require 
just truncating the heavy bosonic states. For Bose-Bose 
mixtures like 87 Rb- 133 Cs (4^|, the modification requires 
allowing for the superfluid order to occur. 



k B T 2 = 



A 



7ry t - 1 [p/(Q) + p6(o)] 



(62) 



Here pf{0) is the density of fermions at the center of the 
trap and pb(0) the density of bosons at the center of the 
trap. 

In Fig. we show T\ and T2 as a function of input 
temperature. Overall, we find very good agreement be- 
tween T\ and T 2 with the input temperature [Fig. [8] (a)]. 
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